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An  efficient  finite  element  (FE)  scheme  to  deal  with  a  class  of  coupled  fluid-solid  problems 
is  presented.  The  main  ingredients  of  such  methodology  are:  an  accurate  Q1/P0  solid 
element  (trilinear  in  velocities  and  constant  piecewise-discontinuous  pressures);  a  large 
deformation  plasticity  model;  an  algorithm  to  deal  with  material  failure,  cracking  propaga¬ 
tion  and  fragment  formation;  and  a  fragment  rigidization  methodology  to  avoid  the  possible 
numerical  instabilities  that  may  produce  pieces  of  material  flying  away  from  the  main  solid 
body.  All  the  mentioned  schemes  have  been  fully  parallelized  and  coupled  using  a  loose- 
embedded  procedure  with  a  well-established  and  validated  computational  fluid  dynamics 
(CFD)  code  (FEFLO).  A  CSD  and  a  CFD/CSD  coupled  case  are  presented  and  analyzed. 

©  2009  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

In  the  analysis  of  complex  structures  subjected  to  blast  loading  and  large  deformations,  Lagrangian  FE  codes  with  explicit 
time  integration  are  a  necessary  tool  [1,2].  This  class  of  problems  requires  hundreds  of  thousands  of  time  steps,  and  in  some 
cases,  millions  of  elements,  which  forces  the  use  of  computationally  simple  methodologies.  Furthermore,  the  material  failure 
and  the  possible  changes  of  topology  due  to  the  crack  advance  has  to  be  implemented  in  an  efficient  manner.  Expensive 
remeshing  procedures  and/or  cracking  schemes  could  make  the  calculation  of  real  problems  an  impossible  task  in  terms 
of  CPU  time.  Moreover,  the  possible  numerical  instabilities  that  may  appear  as  a  result  of  several  fragment  of  elements  flying 
away  has  to  be  taken  into  account.  For  most  cases,  i.e.,  for  bomb  fragmentation  studies,  it  is  not  the  internal  stresses  of  the 
small  fragments  what  needs  to  be  computed  in  an  accurate  manner,  but  their  size  and  shape.  Hence,  a  fast  rigidization  algo¬ 
rithm  that  maintains  the  volume  of  the  fragments  must  be  available.  All  the  mentioned  schemes  should  be  fully  parallelized 
to  obtain  solutions  in  a  realistic  period  of  time. 

The  main  ingredients  of  an  efficient  methodology  to  deal  with  real  coupled  fluid-solid  blast  problems  will  be  described  in 
this  paper.  First,  a  brief  description  of  a  large  deformation  plasticity  constitutive  material  [17]  is  presented.  Such  model  is 
motivated  by  a  well-understood  micromechanical  picture  of  single-crystal  metal  plasticity.  The  resulting  constitutive  model 
relies  on  a  hyperelastic  characterization  of  the  elastic  material  response,  which  avoids  the  drawbacks  of  the  widely  used  hyp- 
oelastic  models,  i.e.,  the  material  isotropy,  the  nonzero  residual  strain  in  a  closed  “elastic”  cycle,  and  in  general,  the  lack  of  a 
stored-energy  function  to  obtain  the  elastic  stress  tensor.  In  addition,  the  deviatoric  and  volumetric  elastic  stresses  are  func¬ 
tion  of  tensors  which  are  invariant  against  rigid  rotations  (tensors  defined  in  the  convective  coordinates).  This  feature  auto¬ 
matically  makes  the  scheme  to  fulfill  the  objectivity  [17]  principle  (no  special  objective  stress  rates  are  used). 
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From  the  numerical  point  of  view,  the  element  kinematic  variables  are  fully  integrated  (eight  integration  points  per  ele¬ 
ment),  which  avoids  spureos  modes  and  its  immediate  consequence:  the  use  of  artificial  hourglass  viscosities  [13].  On  the 
other  hand,  the  pressure  field  is  interpolated  using  only  one  nodal  point  per  element  (piecewise  constant  pressures  per  ele¬ 
ment).  Such  numerical  feature  removes  the  volumetric  locking  and  accurately  approximates  the  material  incompressibility 
in  the  plastic  range  [14].  Also,  a  pressure  smoothing  is  applied  in  each  time  step,  to  avoid  the  appearance  of  check-board  [14] 
pressure  modes,  which  may  be  presented  in  some  cases  using  the  Q1  / PO  element  (this  element  does  not  strictly  fulfill  the 
incompressible  BB  condition)  [14]. 

In  Section  3,  the  material  fracture  model  is  exposed.  For  this  purpose,  the  elements  conforming  the  FE  mesh  are  checked 
using  some  failure  criteria,  i.e.,  the  maximum  effective  plastic  strain,  or  some  damage  variable.  The  most  failed  node  and  face 
of  the  element  are  then  chosen  to  define  the  failure  plane.  The  new  node  is  created  by  duplicating  the  most  failed  one,  and  by 
redefining  the  connectivity  of  the  elements  located  at  one  side  of  the  failure  plane.  Such  methodology  produces  a  cracking 
propagation  algorithm  that  does  not  require  expensive  remeshing  schemes.  This  last  feature  is  highly  desirable  when  dealing 
with  real  problems,  where  millions  of  elements  and  thousands  of  time  steps  may  be  required.  However,  it  may  fails  to  pro¬ 
duce  real  fragmentation  patterns  for  some  cases,  which  motivates  to  implement  semi-empirical  agglomeration  algorithms 

1. e.,  Mott’s  theory  [27,28]. 

A  simple  rigidization  algorithm  based  on  the  number  of  elements  that  a  fragment  of  material  contains  has  been  imple¬ 
mented  to  avoid  numerical  instabilities  and  to  save  CPU  time.  All  the  mentioned  schemes  have  been  fully  parallelized  and 
coupled  using  a  loose-embedded  procedure  with  the  well-established  and  validated  computational  fluid  dynamics  (CFD) 
code,  FEFL098.  Finally,  A  CSD  and  a  CFD/CSD  coupled  cases  are  presented. 

2.  Constitutive  equation  and  numerical  approximation 

As  mentioned  in  Section  1,  the  plasticity  model  used  in  this  work  is  motivated  by  a  well-understood  micromechanical 
picture  of  single-crystal  metal  plasticity  [17].  The  resulting  constitutive  model  relies  on  a  hyperelastic  characterization  of 
the  elastic  material  response,  which  avoids  the  drawbacks  of  the  widely  used  hypoelastic  models,  i.e.,  the  material  isotropy, 
the  nonzero  residual  strain  in  a  closed  “elastic”  cycle,  and  in  general,  the  lack  of  a  stored-energy  function  to  obtain  the  elastic 
stress  tensor.  Moreover,  the  deviatoric  and  volumetric  elastic  stresses  are  function  of  tensors  which  are  invariant  against  ri¬ 
gid  rotations  (tensors  defined  in  the  convective  coordinates).  This  feature  automatically  enforces  objectivity  [17]  on  the 
scheme. 

Also,  it  can  be  demonstrated  that  the  hypoelastic  models  based  on  commonly  used  objective  stress  rates,  like  the  Jau- 
mann-Zaremba  and  Green-Mclnnis-Naghdi  stress  rates,  do  not  result  in  a  straightforward  generalization  of  the  classical  ra¬ 
dial  return  algorithm  when  one  tried  to  develop  a  J2  plasticity  theory  for  large  deformation  cases  (see  [17]).  The  only  models 
which  result  in  such  a  generalization  are  the  ones  based  on  the  Lie  derivate  of  the  Kirchhoff  stress.  However,  many  authors 
copy  as  it  is  the  structure  of  the  radial  return  map  for  small  deformation  plasticity,  to  use  it  with  hypoelastic  models  that  are 
not  based  on  the  mentioned  derivate.  Our  experience  indicates  that  such  algorithms  violates  in  an  unacceptable  manner  the 
isochoric  (incompressible)  character  of  the  plastic  flow  when  non-structured  meshes  are  used:  This  is,  the  elements  tend  to 
have  negative  volume  and  their  edges  tend  to  cross  at  high  explosive  loads  (pressures).  However,  the  use  of  non-structured 
meshes  is  almost  mandatory  for  this  type  of  problems,  not  only  because  they  allow  the  automatic  generation  of  the  compu¬ 
tational  grid  by  using  a  tetrahedra  mesh  generator  (then  each  tetrahedra  is  divided  in  four  hexahedral  elements  see  Fig.  1 ), 
but  also  because  the  unstructured  topology  of  the  mesh  automatically  includes  a  random  behavior  to  the  fracture  algorithm, 
which  makes  the  simulations  more  realistic.  This  point  is  demonstrated  below  in  the  numerical  simulations. 

Ref.  [17]  presents  a  detailed  description  of  the  hyperelastic  plasticity  model  that  has  been  used  in  this  work.  However,  a 
brief  summary  is  shown  below  for  completeness.  The  stored-energy  function  has  the  form: 


w  =  U(Je)  +  W(be) 

a) 

u(f)  =  ^[lye2-  i)-in/ 

(2) 

W(¥)=  lp(tr[F]-3) 

(3) 

which  results  in  the  following  stress-strain  relationships: 

T  =/p/  +  S 

(4) 

P  =  f  c/e2  -  1)/Je 

(5) 

s  =  fidev[be] 

(6) 

Above,  t  is  the  Kirchhoff  stress  tensor,  p  the  mechanical  pressure  and  s  is  the  deviatoric  part  of  the  stress  tensor.  Je  is  the 
determinant  of  the  elastic  part  of  the  deformation  gradient  tensor  Fe,  k  the  material  bulk  modulus  and  p  is  the  shear  mod¬ 
ulus.  Finally,  be  is  the  volume-preserving  part  of  be  (elastic  left  Cauchy-Green  tensor),  given  by: 

bl=f-2/3FeFeT. 


(7) 
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Fig.  1.  From  left  to  right:  finite  element  mesh  of  ^100,000  hexahedrals;  vertical  cut  showing  the  interior  top  part  of  the  weapon. 


The  yield  condition  is  a  classical  Mises-Hubert  type,  but  formulated  in  terms  of  the  Kirchhoff  stress  tensor: 


/M)  ^ll^ll 


[gy  +  Kot]  ^  0 


(8) 


where  gy  denotes  the  flow  stress,  I<  denotes  the  isotropic  hardening  modulus,  and  a  the  hardening  parameter.  ||s||  is  the 
norm  of  s  =  (repeated  index  sum).  It  is  clear  that  non-linear  hardening  laws  (/((a)),  dynamic  increase  factors  for  gy 
and  empirical  based  models  for  the  flow  stress  (Gy)  are  easily  accommodated  in  the  formulation.  As  an  example,  the  widely 
used  Johnson-Cook  [29,30]  model  for  the  flow  stress  in  metals,  which  takes  into  account  the  effect  of  high  strain  rate  hard¬ 
ening  and  temperature  softening,  was  easily  implemented  in  this  plasticity  framework. 

Finally,  the  plasticity  model  is  completed,  as  usual,  with  the  following  associative  flow  rule,  hardening  law,  Kuhn-Tucker 
and  consistency  conditions: 


Me  =  -|ytr[ft>  (9) 


y^O,  /( t,  a)  <  0,  y/(T,  a)  =  0 

y/(T,a)  =0. 


Above,  Lvbe  is  the  Lie  derivative  of  be  and  y  is  the  consistency  parameter.  A  is  the  material  time  derivative  of  A. 


(11) 

(12) 

(13) 
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2  A.  Discretization  of  the  dynamic  equation 


The  dynamic  equilibrium  equation  written  on  the  spatial  configuration  at  time  t(Qt)  is  given  by: 


pit  -  V  •  a  =  F  in 


(14) 


plus  compatible  initial  and  boundary  condition.  In  Eq.  (14)  p  is  the  density,  a  the  Cauchy  stress  tensor,  u  the  displacement 
field,  ii  the  acceleration  and  F  is  the  body  forces.  All  the  terms  are  evaluated  in  the  actual  configuration. 

The  equilibrium  variational  form  at  the  spatial  configuration  is  obtained  by  multiplying  Eq.  (14)  times  a  test  function 
Su  e  V  that  fulfills  the  imposed  boundary  condition  and  continuity  requirements.  After  integrating  by  parts  the  internal 
stress  term  in  the  usual  manner,  the  following  is  obtained: 


(15) 


where  s(-)  is  the  spatial  symmetric  gradient  operator,  s(Su)  is  nothing  but  the  spatial  strain  rate  tensor,  which  is  obtained 
from  the  standard  FE  kinematic  interpolation  (see  Ref.  [14]  for  details),  t  are  the  surface  tractions  which  are  applied  on 
the  boundary  of  Q  (dQ)  in  the  actual  configuration,  and  rN  is  the  Neumann  part  of  dQ. 

It  is  widely  known  that  by  using  standard  FE  functions  and  equal  interpolation  for  the  kinematic  variables  (displacements, 
velocities,  and  accelerations)  and  the  pressure  field  to  discretize  Eq.  (15)  in  space,  a  lack  of  stability  may  appear  if  the  mate¬ 
rial  becomes  almost  incompressible.  This  is  the  case  for  metals  in  the  plastic  range,  where  it  is  standard  to  assume  that  the 
plastic  flow  is  isochoric  (the  material  flows  preserving  its  volume).  The  mentioned  lack  of  stability  (a  pure  numerical  prob¬ 
lem)  is  demonstrated  by  the  well  documented  pressure  check-board  mode  [13,14]  and  locking  effect.  It  is  also  well-estab¬ 
lished  that  the  both  anomalies  may  be  avoided  by  several  numerical  techniques.  Among  them  are:  uniform  reduced 
integration  [14,15],  selective  integration  [14],  mixed  interpolation  [13,14]  (div-stable  elements),  and  numerical  or  physical 
stabilizations  [3-12].  A  mixed  interpolation  has  been  used  in  this  work:  the  kinematic  variables  are  interpolated  using  the 
standard  Q1  FE  functions  (trilinear  8-nodes  hexahedral  elements),  and  the  pressure  is  interpolated  using  the  constant  piece- 
wise-discontinuous  (PO)  FE  function  (constant  pressure  at  elemental  level).  In  addition,  a  filter  of  pressures  (smoothing)  is 
performed  each  time  step  to  enforce  a  convergent  pressure  field  (free  of  check-board  modes).  This  element  does  not  require 
artificial  viscosity  terms  due  to  the  fact  that  it  is  free  of  hourglass  modes  [13]  (the  deviatoric  stress  terms  are  fully 
integrated). 

3.  Fracture  scheme 

Most  of  the  fracture  schemes  may  be  classified  into  two  groups  (or  a  combination  of  them):  smeared  crack  models  and 
discrete  crack  approaches  [18,19].  The  former,  also  called  continuum  approaches,  are  based  on  the  classical  continuum 
(strain  localization,  smeared  cracking)  or  the  enriched  continuum  [16,20,21]  (gradient  enrichment,  Cosserat  continuum, 
non-local  models).  In  such  models,  the  material  softens  (using  some  constitutive  law)  forming  a  band  where  the  crack  (or 
micro-cracks)  grows.  In  general,  the  topology  of  the  mesh  does  not  change  during  the  simulation. 

In  the  discrete  crack  approaches,  the  mesh  topology  changes  to  follow  the  cracks  [18,22].  These  schemes  explicitly  rep¬ 
resent  the  crack  as  a  separation  of  nodes.  When  some  failure  criterion  is  fulfilled,  the  node  is  redefined  as  two  nodes  and  the 
elements  are  allowed  to  separate.  While  this  produces  a  realistic  representation  of  the  opening  crack,  a  coarse  discretization 
in  the  finite  element  model  results  in  misrepresentation  of  the  propagating  crack  tip.  One  way  to  avoid  this  is  by  using  reme¬ 
shing  procedures  in  the  failure  zone  [22]. 

In  the  present  work  a  discrete  crack  approach  has  been  adopted.  The  main  reason  for  this  choice  is  that  in  coupled  blast 
simulations  big  changes  of  topology  are  expected.  A  smeared  approach  would  not  be  capable  of  realistically  simulating  the 
hundreds  of  fragments  flying  away  from  the  main  solid  body,  and  the  possibility  of  impacting  some  other  surfaces.  Reme¬ 
shing  procedures  have  not  been  implemented  yet  in  the  code,  and  it  is  not  mandatory  as  long  as  fine  enough  meshes  are 
utilized.  In  addition,  expensive  remeshing  procedures  could  render  large  scale  cases  impossible  to  simulate. 

The  discrete  fragmentation  scheme  is  implemented  as  follows:  first,  the  adopted  failure  criterion  is  computed  at  the  inte¬ 
rior  of  the  elements.  The  most  failed  element  e  (element  that  violates  the  failure  criterion  with  the  highest  value)  is  then 
chosen.  Among  the  faces  that  are  attached  to  e,  the  one  with  maximum  tensile  stress  will  be  detached.  This  is,  the  one  that 
fulfills  the  following  expression: 


(n i  ■  <re)  ■  n i  ^  (tij  •  <re)  ■  iij ,  Vij  e  and  i  #  j 


(16) 


where  ae  is  the  stress  tensor  in  the  element  e,  n,  the  unit  exterior  normal  on  the  face  i  (exterior  to  the  element  e)  and  is 
the  set  of  faces  attached  to  e. 

To  detach  the  face  with  the  highest  tensile  stress  ,  each  point  ip  e  (4  points  for  hexahedral  elements)  is  duplicated 
in  the  elements  containing  the  topological  perpendicular  edge  to  the  face.  In  other  words,  the  nodal  point  ip  belonging  to  the 
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face  #7  is  duplicated  in  all  the  elements  containing  the  edge  ( ipjp ),  where  jp  is  the  point  that  fulfills  the  following  require¬ 
ments:  (1 )  jp  belongs  to  the  set  of  nodes  of  element  e;  (2)  jp  is  connected  to  point  ip;  (3)  jp  is  not  contained  in  the  set  of  nodes 
of  the  face 

The  algorithm  sketched  above  produces  the  main  failure  scheme  for  weapon  fragmentation  applications:  longitudi¬ 
nal  cracks  consisting  of  very  long  stripes  of  metal  which,  although  theoretically  correct  (tangential  stress  is  higher 
than  the  longitudinal  one),  does  not  yield  the  small  fragment  pattern  observed  in  experimental  tests  [25,26].  For  that 
reason,  a  secondary  fragmentation  algorithm  based  on  a  lot  of  experimental  data  collected  for  more  than  50  years, 
has  been  introduced  to  the  CSD  code.  Such  algorithm  is  based  on  the  Mott’s  theory  of  break-up  of  cylindrical 
“ring-bombs”.  In  [27]  a  detailed  description  of  the  model  can  be  consulted.  Below,  just  a  brief  summary  is  presented 
for  completeness. 

According  to  Mott  [28],  the  average  circumferential  length  of  the  resulting  fragment  is  given  by: 


*o 


r 

V 


(17) 


where  p  and  <jf  are  the  density  and  strength  of  the  bomb  material,  respectively,  r  the  radius  of  the  ring  at  the  moment  of 
fracture,  V  its  radial  velocity,  and  y'  denotes  a  semi-empirical  statistical  constant  determining  the  dynamic  fracture  proper¬ 
ties  of  the  material.  The  value  of  this  last  constant  has  been  calibrated  [27]  using  theoretical  assumptions  and  experimental 
data  to  a  value  of  y'  =  56.  This  has  been  the  value  used  in  this  work  for  all  the  numerical  simulations.  Two  others  important 
relationship  has  been  used  for  relatively  large  fragments.  These  are:  the  representative  values  of  average  aspect  ratios  of 
fragment  lengths  to  circumferential  breadth  l0/x0,  and  the  aspect  ratios  of  circumferential  breadth  to  fragment  thickness 
x0/t0.  The  reported  values  in  the  literature  [27,28]  for  such  relationships  are  between  l0/x0 « x0/t0 « 2.5  and 
lo/Xo  ttXo/to  «  5.0.  In  this  job  a  value  of  4.0  has  been  used.  However,  some  additional  random  effects  are  introduced  to 
the  fragment  sizes.  These  are:  the  unstructured  character  of  the  used  meshes  as  the  reader  can  observe  in  Fig.  1,  and  a  uni¬ 
formly  distributed  random  variable  which  will  be  clarified  later. 

Therefore,  for  weapon  fragmentation  cases  the  final  fracture  algorithm  is  as  follows: 


(1)  Compute  the  element  with  the  highest  value  of  the  failure  criteria  emax  (i.e.,  the  element  with  maximum  effective 
plastic  strain). 

(2)  Obtain  the  circumferential  length  of  the  possible  resulting  fragment  using  expression  (17). 

(3)  Obtain  the  length  of  the  possible  resulting  fragment  using  l0/x :0  =  4. 

(4)  Compute  the  average  failure  criteria /top  in  a  ring  of  elements  of  size  l0  above  emax. 

(5)  Compute  the  average  failure  criteria  fbot  in  a  ring  of  elements  of  size  l0  below  emax. 

(6)  If  /^p  greater  than  the  value  the  material  can  stands,  the  top  ring  will  be  fragmented  using  the  relations 
lo/Xo  =  Xo/to  =  4  *  r,  where  r  is  a  uniformly  distributed  random  number  between  0.625  and  1.25,  that  introduces 
a  random  effect  to  take  into  account,  in  some  way,  the  spatial  variability  of  the  material  strength. 

(7)  If  fbot  greater  than  the  value  the  material  can  stands,  the  bottom  ring  will  be  fragmented  using  the  relations 

Io/Xq  =  Xo/to  =4  *  r. 


The  above  fragmentation  algorithm  produces  results  that  compares  remarkably  with  the  experimental  data. 

3.1.  Rigidization  scheme 

After  the  failure  algorithm  is  performed  and  some  topology  changes  are  detected,  the  possible  detachment  of  fragments  of 
elements  from  the  main  body  is  checked.  The  procedure  is  as  follows:  one  element  e  is  marked  with  a  fragment  identification 
number  ID.  After  that,  the  elements  surrounding  e  are  marked  with  the  same  ID,  the  elements  surrounding  the  elements  that 
surrounds  e  are  also  marked,  and  the  procedure  is  successively  repeated  until  no  more  elements  are  left.  So,  the  fragment 
number  ID  has  been  identified.  Other  unmarked  elements  are  then  chosen  to  attempt  to  identify  other  fragments  by  repeat¬ 
ing  the  whole  procedure.  The  fragment  identification  algorithm  is  finished  when  all  the  mesh  elements  have  been  marked  as 
belonging  to  some  fragment. 

After  that,  all  the  fragments  with  a  number  of  elements  that  is  less  than  a  user  defined  value  are  rigidized.  This  is,  the 
internal  forces  of  the  elements  belonging  to  such  fragments  are  set  to  zero,  and  their  external  forces  are  averaged  and  applied 
as  a  constant  force  to  all  the  fragment  elements.  The  fragment  will  then  behave  as  a  rigid  body. 

Finally,  it  is  important  to  remark  that  the  rigidization  algorithm  is  activated  only  when  the  detailed  stress  distribution 
into  the  small  fragments  is  not  required,  but  their  velocities  and  masses.  This  is  the  case  for  most  of  the  weapon  fragmen¬ 
tation  problems. 


4.  Fluid/Solid  coupling 

The  fluid-structure  interaction  is  simulated  by  coupling  a  CSD  solver  and  a  CFD  solver  in  a  staggered  manner  (see  Box  1 
for  details),  and  by  using  an  embedded  approach  for  the  fluid-solid  interaction  (FSI)  treatment.  In  our  scheme,  the 
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computational  mesh  for  the  solid  body  (for  the  CSD  solver)  is  totally  embedded  into  the  fluid  mesh  (mesh  for  the  CFD  code). 
This  is,  the  fluid  mesh  covers  the  entire  solid  body  (or  bodies),  and  it  is  not  conforming  (not  fitted)  to  the  solid  boundary  faces 
(fluid  mesh  edges  are  intersected  by  solid  boundary  faces). 


Box.  1.  Description  of  the  CFD/CSD  coupling  algorithm. 

Read  meshes;  read  initial  and  boundary  conditions  for  both  solvers; 

CSD  code  sends  the  external  solid  faces  (£f)  and  initial  velocities  (v&>)  on  £?  to  CFD  solver; 

tcFD  <—  0;  tcsD  <—  0; 

DO  WHILE  (tcsD  <  t final  and  tcFD  <  f final) 

CFD  code  computes  the  fluid  elements  intersected  by 

CFD  code  solves  the  fluid  from  t  to  t  +  Stc fd  using  an  embedded  technique; 

CFD  code  sends  the  pressures  P  on  5^  to  CSD  code; 
tcFD  t  4*  <5tcFD  J 
DO  WHILE  (tCSD  <  tern) 

CSD  code  solves  the  structure  (updates  v?)  loading  with  P; 
fcSD  t  +  <5tCSD; 

Update  the  external  faces 
END  DO 

CSD  code  sends  the  new  external  solid  faces  (Sf)  and  velocities 
{y<r)  to  CFD  solver; 

END  DO 


The  staggered  algorithm  and  FSI  treatment  may  be  read  as  follows:  The  simulation  starts  and  all  the  data  for  both  solvers 
are  read.  The  CSD  solver  computes  the  external  faces  mesh  of  the  solid  body,  and  sends  its  coordinates,  connectivities,  and 
velocities  to  the  CFD  solver.  Then,  the  CFD  solver  computes  the  fluid  mesh  edges  that  are  intersected  by  the  solid  faces,  and 
imposes  the  right  flow  boundary  and  initial  (previous  time  step)  conditions  at  the  intersecting  points.  This  is,  the  velocities 
are  set  to  the  solid  faces  velocities  (or  only  to  their  normal  components  for  Euler  type  flows),  and,  in  general,  all  the  variables 
are  interpolated  to  the  mentioned  intersecting  points,  and  to  the  fluid  nodes  which  were  inside  the  solid  body  at  the  previous 
time  step,  if  it  is  required.  This  last  procedure  may  be  performed  by  a  simple  first  order  approach,  or  by  a  second  order 
scheme  which  uses  point  mirroring  and  unknowns  interpolation  with  gradient  reconstruction.  A  detailed  description  of 
the  FSI  (fluid-solid  interaction)  treatment  may  be  consulted  in  Refs.  [25,26]. 

After  the  flow  field  is  solved  taking  into  account  the  blocking  effect  of  the  solid  faces  (boundary  and  initial  conditions 
described  above),  the  flow  pressures  are  interpolated  (applied)  over  such  a  solid  faces.  Then,  the  CSD  solver  updates  the 
structure  kinematic  variables,  internal  stresses,  strains  and  geometry  with  the  new  loads  (new  pressures  over  the  solid 
faces).  Finally,  the  solid  faces  are  updated  and  sent  back  to  the  flow  solver  to  repeat  the  whole  procedure  for  a  new  time  step. 
A  detailed  algorithm  is  shown  in  Box  1. 

The  stability  of  the  CSD  and  CFD  solvers  is  guaranteed  by  choosing  different  time  step  sizes  for  each  code.  The  time  step 
size  for  the  explicit  CSD  code  (<5£csd)  is  computed  in  a  standard  manner,  based  on  the  element  types  and  sizes  [2].  A  standard 
second-order  explicit  Newmark  scheme  (the  widely  called  central  difference  scheme)  is  used  to  integrate  the  dynamic  equi¬ 
librium  equation  in  time  [14].  The  explicit  CFD  code  also  computes  its  own  time  step  size  (<5£Cfd)  using  a  standard  Courant 
criterion  [23]. 

Finally,  it  is  important  to  remark  that  both  the  CFD  and  the  CSD  solvers  are  parallelized  for  shared  memory  architectures 
using  the  Open  MP  libraries.  The  CSD  solver  uses  a  coloring  algorithm  [23]  to  compute  all  the  operations  over  the  elements  in 
a  parallel  manner,  and  a  renumbering  technique  to  avoid  cache  misses.  Most  of  the  cracking  algorithm  is  also  parallelized: 
the  failure  criterion  computation  over  the  elements,  its  smoothing  from  the  element  to  the  nodes,  and  the  computation  of 
imax  (see  Section  3)  are  operations  that  might  be  perfectly  parallelized.  On  the  other  hand,  the  node  duplication  and  element 
reconnection  are  performed  in  a  scalar  way.  This  drawback  does  not  have  much  impact  in  the  total  cost  of  the  simulation 
because  such  operations  involve  a  very  small  group  of  elements.  Therefore,  its  CPU  time  can  be  almost  neglected  in  compar¬ 
ison  to  other  elemental  tasks:  internal  force  computations,  constitutive  law  evaluation,  etc. 


5.  Numerical  examples 

Several  numerical  examples  for  elastic  and  plastic  materials  were  performed  to  validate  the  implemented  solid  element 
[24].  Among  them  were:  the  bar  impacting  on  a  rigid  wall  (typical  benchmark)  and  the  necking  of  a  circular  bar  [17].  All  the 
numerical  results  showed  great  agreement  with  the  ones  reported  by  other  authors. 

Due  to  the  scope  of  this  work,  only  real  coupled  test  cases  are  of  interest.  Hence,  the  fragmentation  of  a  generic  weapon 
inside  a  chamber  with  very  strong  walls  is  presented  to  demonstrated  the  described  methodology. 
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Fig.  2.  Comparison  of  agglomeration  schemes.  From  left  to  right  and  top  to  bottom:  long  strips  and  one-element  fragments  at  100.0  |is  using  the  node 
disconnection  scheme  as  it  is;  same  result  at  200.0  ps  using  the  node  disconnection  scheme;  fragmentation  using  the  Mott’s  approach  at  100.0  pis; 
fragments  at  200.0  ps  using  the  Mott’s  approach. 
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The  geometry  and  some  details  of  the  finite  element  mesh  for  the  weapon  can  be  observed  in  Fig.  1.  The  walls  of  the 
chamber  were  assumed  rigid:  all  its  kinematic  and  dynamic  variables  were  prescribed  to  zero.  The  adopted  weapon  material 
properties  are  presented  as  follows:  Density  p  =  7.76  g/cm3,  Young  modulus  E  =  6.415  x  1012  dy/cm2,  Poisson  ratio  v  =  0.3, 
yield  stress  oy  =  1.283  x  1010  dy/cm2,  isotropic  hardening  modulus  K  =  8.657  x  1010  dy/cm2  and  maximum  effective  plas¬ 
tic  strain  for  failure  e^ax  =  0.3. 


Fig.  3.  Failure  evolution  and  fragment  formation  using  the  node  disconnection  scheme  and  a  structured  mesh.  From  left  to  right  and  top  to  bottom:  fracture 
at  100,  200,  300,  and  400  pis. 
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To  compare  the  straightforward  failure  scheme  (node  disconnection  as  it  is)  and  the  one  used  the  empirical  agglomera¬ 
tion  algorithm  given  by  Mott’s  expression  (17),  a  cheap  stand-alone  simulation  were  performed.  Basically,  the  flow  field 
around  and  inside  the  weapon  was  not  solved,  and  just  a  pressure  front  traveling  inside  the  bomb  from  the  top  to  the  bottom 
were  prescribed.  The  pressure  peak  was  set  to  1011  dy/cm2  and  the  front  velocity  to  7.2  x  105  cm/s.  The  decay  of  the  pres¬ 
sure  front  was  simulated  with  an  exponential  function  of  the  form: 


P  =  PPeake~°U  (18) 

where  d  is  the  distance  from  the  faces  above  the  pressure  front  (the  front  is  traveling  from  the  top  to  the  bottom  of  the 
weapon),  to  the  position  of  the  peak  pressure.  The  pressure  of  the  faces  below  the  peak  pressure  were  set  to  zero. 


Fig.  4.  From  left  to  right  and  top  to  bottom:  CSD  (solid  metal)  velocities  at  12.6  pis,  vertical  cut  of  fluid  velocities  at  12.6  |is,  CSD  velocities  at  37.5  |is,  vertical 
cut  of  fluid  velocities  at  37.5  (is  (CFD  velocity  values  are  hidden  for  clearance  reasons). 
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In  Fig.  2  (top)  the  fragmentation  pattern  using  the  un-modified  failure  scheme  is  shown.  One  can  observe  longitudinal 
long  strips  of  material  plus  some  small  one-element  fragments  (elements  that  has  been  disconnected)  which,  as  it  was 
previously  mentioned,  are  theoretically  correct.  However,  the  final  sizes  and  mass  distribution  of  the  fragments  does  not 
agree  with  the  experience  for  this  type  of  weapons  (see  Fig.  2  top  right). 

Using  the  agglomeration  technique  based  on  the  Mott’s  expression  (17),  the  weapon  fragmentation  seems  to  show  a  bet¬ 
ter  agreement  with  the  observed  experimental  data.  In  Fig.  2  (bottom)  the  evolution  of  the  metal  failure  using  this  scheme  is 
presented. 

For  additional  tests,  a  structured  mesh  was  generated  to  check  the  un-modified  (simple  node  disconnection)  failure 
scheme  for  this  type  of  meshes.  In  Fig.  3  the  fracture  evolution  for  the  structure  grid  is  shown.  Even  though  the  fragment 
formation  improved  with  respect  to  the  un-structured  grid,  the  final  fragments  does  not  look  real.  They  are  too  well  shaped, 
due  to  the  fact  that  the  randomness  introduced  by  the  unstructured  grids  is  lost. 
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Fig.  5.  From  left  to  right  and  top  to  bottom:  CSD  velocities  at  62.5  ps,  vertical  cut  of  fluid  velocities  at  62.5  jis,  CSD  velocities  at  87.5  jis,  vertical  cut  of  fluid 
velocities  at  87.5  ps  (CFD  velocity  values  are  hidden  for  clearance  reasons). 
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Based  on  the  previous  results,  a  coupled  simulation  was  performed  using  the  Mott’s  scheme  for  the  fragmentation  pro¬ 
cess.  The  simulation  was  performed  using  the  algorithm  presented  in  Box  1.  An  explosive  was  ignited  (point  detonated)  in¬ 
side  the  weapon  (at  the  top),  generating  a  pressure  front  that  travels  along  the  case,  and  producing  the  fragmentation 
evolution  showed  in  Figs.  4-8.  The  results  can  be  observed  to  be  in  excellent  qualitative  agreement  with  experimental  data. 
The  fragment  velocities,  sizes  and  shapes  are  very  similar  to  those  reported  by  real  tests.  It  can  also  be  noticed  that  the  top 
part  of  the  weapon  is  detached  from  the  body  at  an  early  fragmentation  stage  (something  commonly  observed  in  experi¬ 
ments).  In  the  same  Figs.  4-8  a  vertical  cut  of  the  fluid  velocities  evolution  may  also  be  observed.  The  effect  of  the  embedded 
fragments  on  the  fluid  velocities  seems  to  be  well  captured,  which  validates  the  used  coupled  approach.  Fig.  9  presents  a  cut 
of  the  fluid  pressures  and  of  the  pressures  over  the  solid  faces  at  different  steps.  This  Figure  demonstrates  the  extrapolation 
of  such  a  variable  from  the  flow  field  to  the  embedded  solid  structure. 
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Fig.  6.  From  left  to  right  and  top  to  bottom:  CSD  velocities  at  112.5  |is,  vertical  cut  of  fluid  velocities  at  112.5  |is,  CSD  velocities  at  137.5  jis,  vertical  cut  of 
fluid  velocities  at  137.5  |is  (CFD  velocity  values  are  hidden  for  clearance  reasons). 
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Fig.  7.  From  left  to  right  and  top  to  bottom:  CSD  velocities  at  162.5  jis,  vertical  cut  of  fluid  velocities  at  162.5  pis,  CSD  velocities  at  187.5  pis,  vertical  cut  of 
fluid  velocities  at  187.5|is  (CFD  velocity  values  are  hidden  for  clearance  reasons). 


Finally,  in  Fig.  10  the  fragment  mass  distribution  is  shown.  The  agreement  is  good  with  experimental  data  for  this  type  of 
weapons.  Fig.  1 1  shows  the  evolution  of  the  fragment  energy.  Again,  it  fits  very  well  the  experience  for  this  type  of  weapons: 
most  of  the  explosive  energy  is  transformed  in  kinetic  fragments  energy. 

As  a  final  remark,  it  can  be  observed  that  all  the  fragments  (including  the  top  and  bottom  parts  of  the  weapon)  are  effec¬ 
tively  rigidized  and  they  fly  alone  free  of  instabilities. 


6.  Conclusions  and  final  commentaries 

An  efficient  CSD/CFD  scheme  for  coupled  blast  simulations  was  presented.  It  is  based  on  a  fully  integrated  solid 
hexahedral  element  with  mixed  interpolation  (Q1/P0:  trilinear  in  velocities  and  constant  piecewise-discontinuous 
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Fig.  8.  From  left  to  right  and  top  to  bottom:  CSD  velocities  at  212.5  |is,  vertical  cut  of  fluid  velocities  at  212.5  |is,  CSD  velocities  at  237.5  |is,  vertical  cut  of 
fluid  velocities  at  237.5  jis  (CFD  velocity  values  are  hidden  for  clearance  reasons). 


pressures),  a  well-founded  large  deformation  plasticity  constitutive  model,  and  in  some  fast  procedures  to  treat  the  mate¬ 
rial’s  fracture.  Moreover,  all  the  algorithms  are  fully  parallelized:  the  CPU  speed  up  obtained  in  the  numerical  examples 
was  quasi-optimum.  One  coupled  numerical  example  was  presented  to  demonstrate  the  overall  scheme,  which  shows  very 
realistic  results. 

Furthermore,  it  was  found  during  the  simulations  that  the  use  of  hypoelastic  schemes  to  treat  the  plasticity  in  metals, 
combined  with  widely  used  objective  stress  rates  (Jaumann-Zaremba,  Green-Mclnnis-Naghdi,  etc.)  and  unstructured 
meshes,  violates  in  a  unacceptable  manner  the  incompressible  character  of  the  plastic  flow.  This  produces  negative  volume 
elements  during  the  simulation,  edge  crossing  and,  therefore,  very  poor  results.  Hence,  it  is  highly  recommended  the  uti¬ 
lization  of  hyperelastic  models  and  convective  formulation  to  treat  the  plastic  material  behavior  in  unstructured  FE 
meshes. 
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Fig.  9.  From  left  to  right  and  top  to  bottom:  Interpolated  pressure  over  CSD  faces  at  87.5  (is  (vertical  clipping  to  show  the  interior  solid  faces  of  the  weapon), 
vertical  cut  of  flow  pressures  at  87.5  |is,  interpolated  pressure  over  CSD  faces  at  187.5  |is,  vertical  cut  of  flow  pressures  at  187.5  ps. 


Finally,  it  is  important  to  remark  that  the  described  CSD  scheme  is  being  used  to  solve  real  coupled  blast  cases.  The 
obtained  results  for  the  large  scale  applications  show  to  be  realistic,  and  the  used  CPU  time  is  also  very  encouraging.  To  give 
an  idea  of  the  CSD  code  speed,  a  simulation  using  a  mesh  of  1  x  105  3D  Ql/PO  CSD  elements,  approximately  spent  0.5  s  of 
CPU  on  1 6  ALTIX  SGI  processors. 
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Fig.  10.  Fragment  mass  distribution. 
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Fig.  11.  Fragment  energy  as  a  fraction  of  the  total  flow  (CFD)  energy. 
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